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Vortex states of weakly-interacting Bose-Einstein condensates confined in three-dimensional rotating 
harmonic traps are investigated numerically at zero temperature. The ground state in the rotating 
frame is obtained by propagating the Gross-Pitaevskii (GP) equation for the condensate in imaginary 
time. The total energies between states with and without a vortex are compared, yielding critical 
rotation frequencies that depend on the anisotropy of the trap and the number of atoms. Vortices 
displaced from the center of nonrotating traps are found to have long lifetimes for sufficiently large 
numbers of atoms. The relationship between vortex stability and bound core states is explored. 
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The recent experimental achievement of Bose-Einstein 
condensation (BEC) in trapped ultracold atomic va- 
pors has provided a unique opportunity to inves- 
tigate the superfluid properties of weakly-interacting di- 
lute Bose gases. Mean-field theories, which are usually 
based on the Bogoliubov approximation or its finite- 
temperature extensions yield an excellent descrip- 
tion of both the static and dynamic properties of the 
confined gases Q. These theories also predict that a 
continuum Bose condensate with repulsive interactions 
should be a superfluid, which can exhibit second sound, 
quantized vortices, and persistent currents. While there 
exists some evidence for second sound in trapped con- 
densates pO[ |, vortices in these systems have never been 
observed despite considerable experimental effort [ pT| . 
Numerous techniques for the generation of vortices have 
been suggested, including stirring the condensate with a 
blue detuned laser adiabatic population trans- 
fer via a Raman transition into an angular momen- 
tum state [ p^ljT^ ], spontaneous vortex formation during 
evaporative cooling ||l^ , |l7|] , and rotation of anisotropic 
traps H]. 

Several studies of vortex stability recently have been 
carried out ||l9|-|2^ . The free energy of a singly-quantized 
vortex attains a local maximum when the vortex is cen- 
tered in a stationary trap jl^. In the presence of dis- 
sipation, such a vortex would migrate to the edge of 
the trap and eventually disappear |20). It has been sug- 
gested 1^,^ that this mechanical instability may be re- 
lated to a bound state in the vortex core, corresponding 
to a negative-energy 'anomalous' dipole mode found nu- 
merically in the vortex state at low densities ||2^. It is 
possible to stabilize singly-quantized vortices by rotating 
the trap at an angular frequency f2. When is larger 
than the 'metastability' frequency fio, infinitesimal dis- 
placements of the vortex no longer decrease the system's 
free energy, and the vortex becomes locally stable; above 



the critical frequency > r^o, the vortex state becomes 
the ground state of the condensate . 

While rotation of the confining potential has been pro- 
posed as a method to both nucleate and stabilize vortices 
in trapped Bose gases, the relevant critical angular fre- 
quencies are not presently known. In the present work, 
numerical results are obtained for Bose-condensed atoms 
confined in three-dimensional rotating anisotropic traps 
at zero temperature. The critical frequency f2c is found 
to increase with the degree of anisotropy in the plane 
of rotation. In order to nucleate vortices, however, the 
trapped gas must be rotated either more rapidly than 
fic, or at temperatures above the BEC transition p6| , p7[ . 
While vortices in nonrotating traps are found to be al- 
ways unstable, for large numbers of atoms their lifetimes 
can be very long compared with a trap period. 

The trapped Bose condensate, comprised of A^o 
repulsively- interacting Rb atoms with mass M = 1.44 x 
10~^^ kg and scattering length a ~ lOOag = 5.29 nm [^ , 
obeys the time-dependent Gross-Pitaevskii (GP) equa- 
tion in the rotating reference frame Pq]: 



V'(r,r), (1) 



where the trap potential is t4 = ^ (a;^ + a^y^ + 
the Hartree term is Vh = 47r?7|-0p, and the condensate 
is rotated about the z-axis at the trap center. The ef- 
fects of gravity (along z) are presumed negligible. The 
constant rotation at frequency induces angular mo- 
mentum per particle given by the expectation value 
of Lz ^ i [ydx — xdy). The trapping frequencies are 



{io.:c,ujy,ujy^) = ujx{\,a,l3) with 



27r X 132 rad/s. 



a > 1, and /3 = -\/8 [|l|. Choosing the condensate 
to be normalized to unity yields the scaling parameter 
?7 — Nqu/cIx. Note that energy, length, and time are 
given throughout in scaled harmonic oscillator units hux, 

2t:/ll}x ~ 7.6 ms, 



dx = ^Jn/Mu, 
respectively. 
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The ground-state of the GP equation is found within 
a discrete- variable representation (DVR) |27| by imagi- 
nary time propagation using an adaptive stepsize Runge- 
Kutta integrator. A total of between 40 000 and 130 000 
DVR points of a Gauss-Hermite quadrature are used, and 
all calculations are performed on a standard worksta- 
tion. The stationary ground state in the rotating frame 
is found by setting t = it and solving the diffusion equa- 
tion: 



a^V(r,f) = -(i/-A')V'(r,f), 



(2) 



where H is the GP operator appearing on the right side 
of Eq. (|l|) and /i is the chemical potential. The conden- 
sate wavefunction is assumed to be even under inversion 
of z, and is initially taken to be the vortex- free Thomas- 
Fermi result, which is the time-independent solution of 
Eq. (H), neglecting the kinetic energy operator and L^. 
A vortex is generated by imposing one quantum of cir- 
culation K, a 27r-winding of the phase around the z-axis, 
on the condensate wavefunction at f = 0. At each imag- 
inary timestep, the chemical potential /i is readjusted in 
order to preserve the norm of the wavefunction (i.e. the 
value of Nq). The propagation continues until the right 
side of Eq. (||) is equal to a tolerance 5 < 10~^° defining 
the error in the dimensionless chemical potential. Sta- 
tionary solutions are verified by subsequently integrating 
in real time; any deviations from self-consistency would 
be made manifest by collective motion. 
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FIG. 1. The condensate density with a singly-quantized 
vortex at the origin (solid lines) , is shown projected along the 
X and y-axes. The vortex-free case is shown as dashed lines 
for comparison. Parameters for both cases are A'^o — 10^ and 
a = 2, yielding ^ x 0.12d^. 

The solution to the GP equation for a vortex located 
at the center of a nonrotating anisotropic trap contain- 
ing Nq = 10^ atoms is shown in Fig. |l|. In general, the 
vortex core is found to become decreasingly anisotropic 
as A^o increases. Furthermore, the condensate density 



preserves its overall vortex-free shape far from the ori- 
gin except for a slight overall bulge in order to preserve 
the norm. The structure of the vortex indicates that the 
healing length ^ is governed largely by the local density 
and is only weakly dependent on trap geometry. In the 
Thomas-Fermi (TF) approximation, which is valid for 
large Nq, the healing length scales with the TF i-axis 
radius R = {15al3r])^/'^d^ as ^ ~ {d^/R)d^ @. Indeed, 
the numerics clearly indicate that the mean vortex core 
radius (approximately dx at low densities) shrinks very 
slowly as the TF limit is approached. 

A superfluid subjected to a torque will remain purely 
irrotational until the critical frequency fJc is reached, at 
which point it becomes globally favorable for the system 
to contain a vortex with a single quantum of circulation 
K. In cylindrically- symmetric systems where the Hamil- 
tonian commutes with L^, the circulation and angular 
momentum (with quantum to) are identical; in the ro- 
tating frame, the free energies of the m ^ 0-states are 
shifted by toSI, and f2c is simply the difference in energy 
between the m = 1 and to = states (divided by K). 
In fully anisotropic traps, however, even the k = state 
is shifted, so the applied in Eq. (||) must be increased 
until the free energy curves cross. It is straightforward to 
extend the TF estimate of to include a small devi- 

ation from cylindrical symmetry [ p^ ; neglecting the shift 
of the vortex- free chemical potential (valid for a ~ 1), 
one obtains: 



(3) 



Fig. g shows the critical frequencies for the global sta- 
bility of a vortex with k = 1 at the trap center. For all 
geometries, the critical frequency drops monotonically as 
A^o is increased. For a given number of atoms, the value 
of increases with in-plane anisotropy, similar to the 
behavior found for liquid helium in rotating elliptical con- 
tainers . The energy of vortex formation must exceed 
that of the irrotational velocity field, which is finite for 
a vortex-free condensate in a rotating anisotropic trap. 
The TF result (||) agrees well with the numerical data in 
its regime of validity a ^ 1, though it tends to slightly 
overestimate the value of i7c- 

While ilc provides the criterion for the global stability 
of a vortex, it does not necessarily indicate the critical 
frequency for vortex nucleation. When initially vortex- 
free condensates are placed in anisotropic traps rotating 
at a frequency SI < cj^;, the velocity field of the station- 
ary solution is found to be irrotational even for 3> rJc- 
In a harmonic trap with smooth edges, it is not clear if 
there exists any suitable locus for vortex formation. The 
vortices are most likely to originate at the condensate 
surfaces normal to the axis of weak confinement, where 
the local critical velocity is small but the tangential 
superfiuid velocity in the laboratory frame is largest p9[ . 
While these issues are beyond the scope of the present 



2 



issue of vortex stability, there is evidence that multiple 
vortices appear at higher frequencies For smaller 

A^Oi it would likely be easier to generate a vortex ex- 
perimentally by rotating the anisotropic trap before the 
condensate is cooled below the BEC transition [n6 17|. 
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FIG. 2. The critical frequency for the global stabilization of 
a singly-quantized vortex at the origin is given as a function 
of the number of atoms in the condensate A'^o and in-plane 
anisotropy a. Circles, squares, and triangles show numeri- 
cal results for a = 1,2,3; bold solid, short dashed, and long 
dashed lines are the corresponding TF estimates, respectively. 
The thin solid line represents the frequency of the anomalous 
mode (shown as a positive value) in the vortex state for a = \. 

When a > 1, the angular momentum per particle /k is 
a nontrivial function of A'o, a, and fi. In a nonrotating 
system with unit vorticity, increases with Nq. In the 
absence of a vortex, 1^ is finite for a given f2, and in- 
creases with a; the superfluid velocity Vs can be locally 
appreciable but still remain irrotational V x Vs = 0. At 
the critical frequency, the difference between l\ and is 
always less than unity; for the most extreme case consid- 
ered here, a system with A^o = 10^ and q = 3 rotating 
at Sic — 0.1407^:, one obtains l\ — 2.63?i and — 1.77?i. 
As a ^ cxD, the angular momentum approaches that of 
a non-superfluid TF cloud ~ Ish^ with 'solid-body' 
moment of inertia /gb = jMR^. 

An anisotropic harmonic oscillator potential becomes 
unconfining when it is rotated at a frequency between 
the smallest and largest trapping frequencies. Since flc 
exceeds lUx for sufficiently large a, there exists a critical 
minimum number of condensed atoms Nc able to sup- 
port a vortex. The value of Nc increases with a and is 
given by the intercept of the Vl/ujx = 1 hne in Fig. ^. In 
cylindrically-symmetric systems Nc — 1, since in the ro- 
tating frame the free energies for all the m-states become 
degenerate dX = lOx- In the limit of extreme anisotropy 
a — > oo vortices can never be stabilized. 

It should be noted that states with vortices at the cen- 
ter of anisotropic harmonic traps are found to be station- 



ary solutions of Eq. (g) for all values of Nq and f2 > 
considered; such configurations do not appear to decay 
in either real or imaginary time. A vortex at the center 
of a nonrotating cylindrical trap increases the system's 
free energy, but is stationary because both the vorticity 
and angular momentum commute with the Hamiltonian; 
in principle, the angular momentum can be eliminated, 
and the free energy reduced, only if this symmetry is 
broken by displacing the vortex from the center. Since 
angular momentum is not conserved in anisotropic traps, 
the apparent vortex stability is likely due to the free en- 
ergy maximum at the trap center |T^. In the absence of 
an external pinning mechanism, any such configuration 
should be unstable against infinitesimal displacements. 

In order to further explore the issue of vortex stabil- 
ity in nonrotating traps, the initial condensate phase is 
wound by 27r a small distance « 0.2^2; from the ori- 
gin of a trap with a = \. For all values of A'o < 10^, the 
condensate wavcfunction rapidly (by f ~ T) converges to 
a metastable solution with a vortex, where the fluctua- 
tions in n become smaller than 5 w 10^^ per timestep 
Af ~ 10"'^T. This wavcfunction subsequently decays 
to the true ground state, but both the real and imag- 
inary time required to do so is found to increase with 
Ao |3|]. To an excellent approximation, the total time 
diverges like f oc A'g^^T; for A^o ~ 10^, the time re- 
quired ~ 30T becomes computationally inaccessible and 
the vortex state becomes numerically indistinguishable 
from stationary. The numerics suggest that while vor- 
tices in nonrotating traps are always unstable against 
off-center displacements, they may be very long-lived. 

The observed j;o > instability of the vortex state 
is likely due to the existence of an 'anomalous' collec- 
tive mode Wa at low densities |jl^,^,^. This dipole 
mode, which has positive norm but negative energy (or 
vice versa) , is associated with a zero angular momentum 
bound state in the vortex core |2^ ; its value corresponds 
to the precession frequency of the vortex relative to the 
cloud |]l8|. Previous numerical calculations Q found 
It^al > for aU No < lO'*. As the core radius shrinks 
with larger Aq, however, the anomalous energy might be 
pushed to zero, yielding long-lived or even stable vortices 
in the TF limit. 

The low-lying excitation frequencies of a nonrotating 
condensate in the vortex state are calculated using the 
Bogoliubov equations For completely anisotropic 

geometries, however, the Bogoliubov operator is too large 
to diagonalize explicitly. Calculations are therefore re- 
stricted to the cylindrical case a = 1, where the vortex 
condensate is = i^i{p,z)e'"^ and the quasiparticle am- 
plitudes u and V are labeled by to, the projection of the 
angular momentum operator L^. The Bogoliubov equa- 
tions are then 
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where 6 ^ -^^f^pf^ - 
O + ^'■^I"'-' , and ui = 



I 

2 ( 



when ei 



2Fh, O' = 
: 0. In the 



p-direction, the points of the DVR grid correspond to 
those of Gauss-Laguerre quadrature, and the kinetic en- 
ergy matrix elements are obtained using the prescription 
of Baye and Heenen . 

The anomalous mode uja, which is labeled by m = 2, is 
shown as a function of TVo in Fig. 0. The results indicate 
that < \uja\ < i^c for all iVg < 10^ considered. Our cal- 
culations suggest that for large numbers of atoms, oja co- 
incides with the metastability rotation frequency f2o dis- 
cussed above; the numerical value oiiiOa is consistent with 
the TF result fio = ff^c [@- Indeed, in the frame of a 
condensate rotating at = w^, the frequency of the vor- 
tex oscillation would be Doppler shifted to zero. Alter- 
natively, it can be shown in both the weakly-interacting 
and TF limits that f^o is also the frequency at which 
the chemical potentials /i for the vortex and vortex-free 
states become equal; in the TF limit, uja vanishes when 
the vacua (or energy zero) for quasiparticle excitations 
for both states coincide. 

Since the anomalous mode corresponds to the preces- 
sion of the vortex about the trap origin, one may make a 
crude estimate of the vortex lifetime r. In the presence 
of dissipation, the vortex will spiral out of the condensate 
after a few orbit periods OJ^^. Assuming that oja = f ^^c, 
then with Eq. (|) one obtains r - N^'^T in the TF limit 
neglecting logarithmic factors. This result is consistent 
with the imaginary time f required to yield the vortex- 
free ground state in the fully three-dimensional numerical 
calculations discussed above. Similar decay times have 
been obtained for solitons and vortices in the presence of 
a small noncondensate component p^ . 

In summary, we have obtained numerically the crit- 
ical frequencies fic for the stabilization of a vortex at 
the center of a rotating anisotropically-trapped Bose con- 
densate. Since fie increases with the in-plane anisotropy 
a — LOy/ujx and the condensate becomes unconfined for 
f2 > Wa;, there is a minimum number of atoms able to 
support a vortex state. Vortices in nonrotating traps are 
found to be unstable against small off-center displace- 
ments, but their decay time diverges with the total num- 
ber of atoms. 
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